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Abstract 

A recent paper [V. L. Campo et ai, Phys. Rev. Lett. 99, 240403 (2007)] has proposed a two-parameter scaling method to 
determine the phase diagram of the fermionic Hubbard model from optical lattice experiments. Motivated by this proposal, we 
investigate in more detail the behavior of the ground-state energy per site as a function of trap size (L) and confining potential 
iyix) = t(x/L) a ) in the one- dimensional case. Using the BALDA-DFT method, we find signatures of critical behavior as a — + oo. 
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The Hubbard model was proposed in the 1960s as a sim- 
plified model of a system of correlated fermions, capturing 
the competition between localisation due to strong inter- 
particle repulsion and itinerant behaviour due to intersite 
hopping. It remains a model of central importance in con- 
densed matter physics; but despite this fact, and decades 
of investigative effort, the model in dimensions d > 1 re- 
mains unsolved, and its phase diagram is still a subject of 
controversy 

Recently, a new line of attack has been proposed pQ : to 
determine the phase diagram of the Hubbard model not us- 
ing theory but using cold-atom experiments which — un- 
like traditional condensed matter systems — should be es- 
sentially perfectly described by the Hubbard Hamiltonian. 
This is for two reasons. Firstly, the lattice in the cold-atom 
case is provided by a laser standing wave (or 'optical lat- 
tice'), which can be made almost perfectly periodic and 
sinusoidal. Secondly, the interaction between the neutral 
atoms should be well described by a contact interaction of 
the Hubbard form [2|3j . Other advantages of such exper- 
iments include the following: the inter-atomic interaction 
strength and hopping amplitude in these systems can be 



tuned over a wide range; the achievable particle numbers 
are much greater than those available in current computer- 
based simulations of fermionic systems; and the total en- 
ergy (for example) is an easily measurable quantity. 

However, a major difference between the solid-state and 
cold-atom realisations of the Hubbard model is the nature 
of the particles' confinement. In the solid-state case, it is of 
the 'hard- wall' type, where the system is homogeneous ex- 
cept for a very high potential step at its edges; in the cold- 
atom case, the confining potential arises from an external 
electromagnetic field, and is almost always harmonic, leav- 
ing no residue of translational invariance. How to extract 
the phase diagram of the homogeneous model from these 
inhomogeneous systems is discussed in Ref. [4], where a 
two-parameter scaling method is proposed, with the two 
parameters L and a representing respectively the size and 
shape of the confining potential. The method involves con- 
sidering different such shapes (i.e. different values of a), and 
determining for each of them a thermodynamic (L — > oo) 
limit of some intensive quantity such as the energy per site. 
Then one analyses the dependence of that quantity on the 
shape of the confining potential and extrapolates (in a) to 
the homogeneous ('hard-wall') limit. 

It has been empirically observed |4] , at least in the one- 
dimensional case, that the rate of convergence to the ther- 
modynamic limit is dependent on the shape of the confin- 
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Fig. 1. Optical lattice and different confining potentials. The optical 
lattice with a high amplitude defines the sites and makes it possible to 
describe the system by means of a Hubbard model, which is however 
inhomogeneous due to the confining potential (thick lines). As the 
exponent a in <(2j tends to oo, the confining potential approaches 
the square well shape with width 2L. 

ing potential. Here, we address that question in more de- 
tail. In the course of so doing, we find that the behavior 
during the second extrapolation, i.e. the extrapolation in 
the shape of the potential towards the homogeneous limit, 
displays signatures of critical behavior. 

2. Model, and previously published results 

We consider the following inhomogeneous one-dimensional 
Hubbard model: 



6 = -*E(4,*W + H.c.) + U^, 



3,a 



(1) 



where Cj a creates a fermion with z-component of spin a at 
site j, t measures the hopping integral between neighbour- 



ing sites, U is the strength of the local interaction between 
two fermions of opposite spins, hj t(7 = cj a Cj ia is the num- 
ber operator for spin a at site j, and V(x) is the external 
potential rendering the model inhomogeneous. The posi- 
tion of the site j is Xj = ja, where a is the lattice parameter 
of the chain, and j takes any integer value. The external 
potential assumes the general form 



V(x) = t 



(2) 



with a > 0, and therefore confines the fermions to a re- 
gion of size ~ 2L, or more precisely, makes the probability 
of finding a fermion at x with \x\ 3> L negligible. Accord- 
ingly, we refer to 2L as the system size, and 2L/a as the 
number of sites, irrespective of the value of a. This allows 
us to define an intensive quantity as the corresponding ex- 
tensive quantity divided by 2L/a. In particular, we will be 



Fig. 2. Dependence of total energy on system size for different trap 
exponents, a = 2 (+), 4 (*), 6 (□), 8 (■) and oo (x). The filling 
is fixed at / = 1. We show results for non-interacting fermions 
(U = 0) on the left and for interacting fermions (U = At) on the 
right. Notice that for finite values of a the thermodynamic limit is 
attained with smaller system sizes than in the a = oo case. Insets: 
the same quantity in the thermodynamic limit plotted against a -1 . 
Figure reproduced from Ref. [4]. 



concerned here with the energy per site, Ea/2L, and with 
the filling / = Na/2L, where N is the number of fermions. 

The above Hubbard model, and its higher dimensional 
versions, can be realised with trapped ultra-cold fermionic 
atoms in the presence of an optical lattice with suitable 
laser intensity and wavelength (see Fig. [I}. The atoms are 
confined by an additional external potential whose expo- 
nent a is usually equal to 2. In optical traps with Gaussian 
laser beams, a suitable superposition of two laser beams 
could eliminate the harmonic part of the potential, leav- 
ing a quartic one. Such trap has already been reported [5]. 
Similarly, superposing three or more laser beams, one could 
generate trap potentials with exponents of 6, 8, and so on. 
Despite the increasing experimental difficulties involved in 
making traps with higher exponents, the two-parameter 
scaling method relies on this possibility. For applications to 
solid-state systems, one would be interested in the thermo- 
dynamic limit [L — > oo) of the above model with a = oo 
(a square well potential of width 2L). 

As Fig. [2] illustrates, given a finite exponent in the con- 
fining potential, simply increasing the size of the trap, L, 
is not enough to recover the homogeneous results. Traps 
with different exponents a have different thermodynamic 
limits of energy per site, and the same must happen with 
other properties. The unavoidable confining potential with 
a finite exponent poses, therefore, an intrinsic obstacle to 
the extraction of the homogeneous model's phase diagram 
from such experiments. To overcome it, one must make a 
sequence of experiments [4] with different exponents and 
then extrapolate the results to get the limit a — > oo. 





Fig. 3. Dependence of Ae = E/2L — too on the system size for 
(a) non-interacting fermions (U = 0) and (b) interacting fermions 
(U = 2t). In both cases, numerical data from systems with filling 
/ = 0.5 were used. For each value of the trap exponent, a, Eq. J3} fits 
the data very well even for quite small systems. Similar behavior is 
found for fillings / = 1.0 and / = 1.5 and also interaction strengths 
U = 4t and U = 8t. 

3. Method, and new results 

In order to achieve a better understanding of the role 
the trap exponent a plays in the dependence of the en- 
ergy per site on the system size, we investigate both non- 
interacting (U — 0) and interacting (U ^ 0) systems. While 
the former require only the diagonalization of a tridiago- 
nal matrix, the latter require methods which can deal with 
the strong correlations coming from the interactions. We 
choose an approach based on Density-Functional Theory 
(DFT) [6] , which has been specially developed to treat the 
one-dimensional inhomogeneous Hubbard model |7|8j . Its 
energy functional is constructed by means of the local den- 
sity approximation based on the exact Bethe Ansatz so- 
lution for the homogeneous model: this approximation is 
termed BALDA. Comparisons between BALDA and Quan- 
tum Monte Carlo (QMC) results demonstrate that BALDA 
gives ground-state energies with an accuracy of a few per- 
cent [5]. 

Our numerical simulations show that as the system size, 
2L, goes to infinity, the system's ground-state energy per 
site approaches its thermodynamic limit according to 



e(L) 



2L 



E(L) 



(3) 



where eoo, £, and 7 all depend on the trap exponent a in a 
way which resembles critical behaviour, with the hard-wall 
case (a = oo) corresponding to the critical point. 

In Fig. [31 we show that the energy per site follows Eq. 
for both non-interacting and interacting systems. Through 
the fitting we can extract the dependence of its three pa- 
rameters on the trap exponent (Fig. Although Fig. [3] 
only displays results for filling / = 0.5, we have found sim- 
ilar behavior for different fillings (/ = 1 and / = 1.5) and 



Fig. 4. Dependence of the thermodynamic limit of energy per site, 
eoo, of the finite-size scaling exponent, 7, and of the correlation 
length, £, on the trap exponent, a, for non-interacting ((a), (b), and 
(c)) and interacting ((d), (e), and (f)) fermions. In the non-interacting 
case, we show results for fillings / = 0.5 (x) and / = 1.0 (o); in the 
interacting (U = 2t) case, we show in addition results for / = 1.5 

(In- 



different interaction strengths (U = At and U = 8t). 

Fig. |4] displays the dependence of the thermodynamic 
limit of energy per site 6^, the finite-size scaling expo- 
nent 7, and the length £ on the trap exponent a. These 
graphs are reminiscent of critical behaviour, with a -1 = 
(i.e. the square-well trap) as the critical point. While the 
thcrmodynamic-limit energy is a continuous function of 
a -1 , the exponent 7 seems to be discontinuous at a" 1 = 0, 
and the length £ seems to diverge as a -1 — > — lending 
support to our use of the term 'correlation length'. 

Note that, in the interacting case, the finite-size scaling 
exponent 7 is a-dependent, and furthermore seems to tend 
to the non-interacting value 7 = 2 as a -1 — ► 0. Also, in 
all finite-a cases we have 7 > 1. This accords with our 
earlier observation that the thermodynamic limit is more 
easily achieved for finite trap exponents. The reason for 
this is likely to be that for finite a the wave functions 
can extend beyond \x\ = L. Thus the system does not 
end abruptly, which in turn avoids strongly perturbing the 
fermionic fluid. 

The apparent discontinuity in 7(a) at a -1 = 0, however, 
remains surprising. In the non-interacting case, we have 
performed calculations for a as large as 200. Looking care- 
fully at Fig. [3K, one can see that in the case a = 200, the 
points corresponding to smaller systems seem to follow a 
line whose 7 is equal to 1; but as the system size becomes 
larger, the energy starts to follow the line with 7 = 2. We 
can, therefore, identify a transition length, Lt, and we ex- 
pect Lt to become larger and larger as the trap exponent 
is increased, diverging when a -1 — > 0. This implies a rela- 
tion between Lt and the correlation length £, though ad- 
ditional computations with larger systems are necessary to 
determine its nature. 

The dependence on filling fraction, /, is also instruc- 
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tive. The finite-size scaling exponent 7 and the correla- 
tion length £ do depend on /; however, the plots indicate 
that this dependence disappears as we approach the criti- 
cal point a' 1 — 0. This is exactly what we would expect 
in a critical scenario: the critical point would be described 
only by its universality class, which would be independent 
of a parameter like the filling fraction. In particular, the 
slopes in Figs. [4J: and [4f are essentially the same for the 
different fillings and give us directly the critical exponent 
v associated to the correlation length. We obtain v ~ 0.51 
for non-interacting fermions and v ~ 0.70 for interacting 
fermions. We must emphasise, however, that this latter is 
a rough estimate only. Computations with larger systems 
will be necessary to determine the true behavior of inter- 
acting fermions near the critical point, even if the BALD A 
method is reliable in such extreme cases. 

4. Conclusions 

In conclusion, we have considered the inhomogeneous 
one-dimensional Hubbard model, studying the evolution 
of the ground-state energy per site towards its thermody- 
namic limit as the system size is increased. Looking at the 
role played by the trap exponent, we have found an appar- 
ently critical behavior, which still needs more detailed char- 
acterization. While the results for non-interacting fermions 
are exact, the results for interacting fermions were obtained 
within the BALDA approximation. 

Since the observed critical behavior has a geometrical 
origin, we expect that such behavior must happen inde- 
pendently of the interaction strength. Our treatment us- 
ing BALDA can be interpreted as a mean field calcula- 
tion, which is able to capture the critical behavior but gives 
inaccurate results for critical exponents, for example. To 
overcome this methodological limitation, one should adopt 
more accurate computational methods, such as QMC [9] 
or Density Matrix Renormalization Group (DMRG) [TO] , 
These are, however, much more computationally intensive 
than a DFT calculation. 

From the analytical point of view, even the non- 
interacting case deserves further exploration. In particular, 
in the absence of interactions it should be easier to find a 
renormalization group approach with a convenient deci- 
mation scheme to extract the critical behavior accurately. 
Once such scheme has been found, given the similarity 
between the interacting and non-interacting results, such 
an approach could be tentatively generalized to the inter- 
acting case. 
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